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Abstract 

The air flow through a propeller-type wind tur- 
bine rotor is characterized by. three-dimensional 
rotating cascade effects about the inner portions 
of the rotor blades and compressibility effects 
about the tip regions of the blades. In the case 
of large rotor diameter and/or increased rotor angu- 
lar speed, the existence of small supersonic zones 
terminated by weak shocks is possible. An exact 
nonlinear mathematical model (called a steady Full 
Potential Equation - FPE) that accounts for the 
above phenomena has been rederived. An artificially 
time dependent version of FPE was iteratively 
solved by a finite volume technique involving an 
artificial viscosity and a three- level consecutive 
mesh refinement. The exact boundary conditions 
were applied by generating a boundary conforming 
periodic computation mesh. 
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speed of sound 

correction to the absolute velocity 
potential at point j 

Jacobian - det (j) 
specific static enthalpy 
rothalpy = h + ^7/^ - 
geometric transformation matrix 
relative Mach number - IV |/a 

= (V r -V r - fiV)/2 

the relative air speed ~ 1 

residue of the full potential equation 
at point j 

position vector in the rotor plane (y,z) 
specific entropy 

orthogonal relative streamline aligned 
coordinate system 

absolute temperature of the air 

artificial time 
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modified contravariant relative velocity 
components in X,Y*2 space 

the relative velocity vector 

the absolute velocity vector - 

ySxr 

global boundary conforming coordinate 
system in computational space 

local coordinates in each computational 
cell in the computational space 
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cp absolute velocity vector potential 

angular rotor speed = const, 
a) over-relaxation factor 


Introduction 

The overall efficiency of the propeller- type 
wind turbines (Fig. 1) presently used for the pro- 
duction of electric energy can be significantly 
improved through various aerodynamic modifications. 
One of the basic disadvantages of all the theories 
used for the design and analysis of blade shapes 
is that they include linear and often semi- 
empirical approxima tions and corrections when ac- 
counting for the nonlinear effects of compressible 
three-dimensional rotating cascade flow. 

The purpose of this paper is to present a 
numerical method for solving an exact three dimen- 
sional full potential equation that models the in- 
viscid, irrotational, homentropic flow of a com- 
pressible fluid through an arbitrarily shaped iso- 
lated rotor. This work is based on the principles 
used in external transonic aerodynamics^ 2 , 3 an d 
represents an extension of the authors research**^, 
in the field of potential transonic axial turbo- 
machinery flows. 


The Mathematical Model 


The derivation of the governing equations is 
based on the following assumptions . Assuming that 
the rain drops, snow flurries, atmospheric ice 
particles, industrial pollutants, sand and dust are 
uniformly distributed throughout the air volume and 
that their total volume and mass are negligible, 
the atmospheric air can be treated as a homocomposi 
tiona^l fluid. In addition, it is necessary to as- 
sume that the air is ipviscid and that the atmo- 
spheric turbulence and disturbances due to the pres 
ence of 'a tower and a ground are negligible. 

Then the full potential equation can be ob- 
tained from the following analysis. If the oncom- 
ing airstream is a xi symmetric with respecc to the 
axis of rotation and the rotor angular speed U is 
constant, the problem becomes a steady one when ex- 
pressed in terms of coordinates (x,y,z) fixed for 
the blade (Fig. 1). 

If V r is the relative velocity vector of the 
air with respect to the blade, r is £he position 
vector in the plane of rotation and {} is the 
angular velocity vector of the oncoming flow, then 

5 = ? r + oy? (1) 

is the absolute velocity vector. Then the sum of 
inertia, centripetal, Coriolis and pressure forces 
is 
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where T is the absolute temperature , S is the 
entropy and rothalpy I is defined^ as 

I = h + ~ (v r -V r - n 2 r 2 ) (3) 

where h is the static enthalpy. In order to be 
able to use a single variable, the so-called abso- 
lute velocity potential function <p(x,y,z) where, 

V = 7q> (4) 

the condition of irrotationality 


and M r *= q r /a the local relative Mach number. 
Type dependent rotated finite differencing 1 >4 was 
used for the discretization of Eq. (11). 

The solution of this steady state equation 
(Eq. (11)) can be obtained as an asymptotic solution 
to an unsteady equation for the large time. 1 ^ In 
order to accelerate the iterative solution process , 
but to avoid small time steps dictated by the numer- 
ical stability condition m the case of a truly un- 
steady FPE, a more general artificially time depen- 
dent form of Eq . (11) was used 
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must be satisfied throughout the flovfield. 

For Eq. (5) to be satisfied, the rothalpy must 
be constant 


71 = 0 (6) 

and the flow homentropic 

7S = 0 (7) 

simultaneously everywhere in the flowfield. 
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The consecutive iteration sweeps were considered as 
steps in an artificial time direction. The artifi- 
cial time dependent derivatives in Eq. (13) were ob- 
tained by a careful arrangement of the absolute 
velocity potentials obtained from the two consecu- 
tive iteration sweeps. For example , 1 >4 the central 
difference approximation of the second derivative 
in the s-direction, evaluated at the point (i,j,k), 
is 


Besides already mentioned assumptions and re- 
strictions, Eqs. (6) and (7) imply that there should 
be no heat transfer between the blades and the air, 
boundary layer should not separate and all possible 
shock waves should be weak. 

The flow at upstream infinity was assumed to 
be uniform, although, according to the analysis 
which was done m the earlier papers,^ it is pos- 
sible to introduce a two-dimensional potential vor- 
tex at x = 
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where 9^ , ^ is the old value of the potential 
evaluated’ during the last iteration sweep, 
is the new value of potential function evaluated 
during the present sweep and (jp^j ^ is a temporary 
value on the line along which the relaxation is ap- 
plied. The last term is defined 1 as 


The continuity equation,' 1 
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■where a is the local speed of sound and, 

Q - \ (V\ - fi 2 r 2 ) ( 9 ) 

can consequently be written^*® in its full potential 
form 

a 2 Vq> - (W-V) \ (VP- 1 ?®) + 2(7p.?)((fi x r).$P) 

- «fi x r)-V)((fi X r)-W) = 0 (10) 

This second order quasilinear partial differ- 
ential equation of the mixed type was numerically 
solved using an iterative successive line over- 
relaxation technique. In order to account for the 
proper domain of influence in the case'' of locally 
supersonic flow, the FPE should be written in its 
canonical form^ 

(4 - - »„) 4 0 

Here, (s,m,n) is an orthogonal coordinate system 
locally aligned with the relative velocity vector 


9. . . = 1 (ot . - <P? . <15) 

I,j,k 03 V I,j,k 1 , J , k / 1,3, k 

where co is the over-relaxation factor. Then 
Eq. (14) becomes 

J.k) 

+ ( V i-l»j,k * ^i-1 j j jk) 
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By adding and subtracting (At)(cp t ). from 

* 1, J j K. 

Eq. (17), we finally^ get 
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All the second derivatives in the full potential 
equation were evaluated using the type dependent 
finite difference approximations. This means that 
it) a . locally subsonic region (M r < 1 ) central dif- 
ferencing (designated by the superscript E) was 
used, while in the case of a locally supersonic 
relative flow (M r > 1) the rotated* upstream differ- 
encing was used (designated by the superscript H). 


D - det 1 Jj (22) 

then the modified contravariant components* of the 
relative velocity vector* are (see Eqs. (1) and (4)) 


V l » dEjT 1 J V 1= dIj]' 1 < %iz i + DiAiJ CP, 
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Computational Mesh 

In order to apply the exact boundary conditions 
(with no approximation on the surface of the blade 
and a rotor hub) it is necessary to generate a com- 
putational mesh that will conform with these irregu- 
lar solid boundary shapes. At the same time this 
mesh should be (preferably) periodic in the 6- 
direction, thus providing for an easy application 
of the periodicity conditions along the arbitrarily 
shaped periodic boundaries (lower and upper bound- 
aries on Fig. 2). 


The full potential equation, which can be written as 

(■ 2 - - 2 ) • *?„) ♦ (AV - ,V ,) ■ . 


after division with (-a 2 ) becomes 
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Figure 2 represents such an irregularly shaped, 
non -orthogonal boundary fitted mesh on one of the 
x, 6; r = const, cylindrical computational planes 
that intersect the blade (see Fig, 3). 

The mesh of Fig. 2 was generated using conform- 
al mapping, elliptic polar coordinates and coordi- 
nate stretchings and shearings * as shown on 
Fig. 4. 

For the purpose of finite differencing, each 
distorted mesh cell is separately mapped into a unit 
cube (Fig. 5) using local isoparametric trilinear 
mapping functions of the general form 

8 

ls iIv i + V (U V (1 + ' 2 V (19) 


< 27) 

The full potential equation can also be written in 
its scalar form* (see Eqs. (8) and (23)) as 

4 V r >Y + W r j2 )’ (L ’r Q ,X + V.Y * W r Q ,Z> = 0 


The evaluation of all the finite difference 
approximations was done by the use of the following 
expressions : 

( b > X \=l/2, Js k = ± 2 ‘ 


where subscript p refers to the value at the 
cube's comer, that is 


X - ±1 
P 


±1 Z = ±1 
P 


and b stands for any of the following: x,y,z,<p. 


Computational Space 


Besides transforming the geometric parameters 
from the physical (x,y,z) into the computational 
(X,Y,Z) space, the same was done numerically with 
the governing full potential equation. 


(\x) i)j±1/2jk 8 ( b i+I,j,k " b i-l,J,k‘ rb i+l,j=l ) k 

* b i-l, j=l,k) 

( b,X ^i 5 j k±l/2 " ® " b i*l j j >k + b i+l, j,k=l 




because the neighboring mesh points are spaced in 
each new direction (X,Y,Z) one unit apart. Here b 
stands for any of the following: x,y,z, Ur* v r> 

W r , Q, <P- 

The analogous formulas are valid for the dif- 
ferences in the Y and Z directions. 

On each column j - MAXY, discretization of 
the FPE leads to the set of nonhomogene ous , non- 


and 
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linear, algebraic equations of the general matrix 
form 


[B] {Cj ) + {^3 = 0 (30) 

Equation (30) was directly solved for the vector of 
corrections (Cj) to the potential <p, where 


C 
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Equation (25) was used for evaluation of tri- 
diagonal coefficient matrix !bJ, while Eq. (28) was 
used for the evaluation of the residual vector (Rj) 
which also incorporates an explicitly added artifi- 
cial viscosity in conservative f ora. ^ >3,4 


programming Considerations 

A computer program called WIND was developed 
on the basis of the previous analysis. The mesh 
generating portion of the code uses about 50% of 
the total high speed memory required by WIND and at 
the same time consumes less than 57c of the total 
CPU time required by WIND* Therefore, in order to 
save on computer storage and at the same time to 
provide the means of separately analysing the geome- 
try, WIND was devided into two separate programs. 

The first portion of this code (WIND-01) generates 
the three-dimensional body fitted computational 
mesh. Actually, the first program generates three 
consecutively refined meshes and stores them on 
separate disks. The second portion of the program 
(WIND-02) reads these (x,y,z) coordinates into the 
high speed memory in such a way that only the data 
from three neighboring cylindrical computational 
planes (see Fig. 3) have to be scored in these 
arrays at one time. 


The hub is defined as a doubly infinite circu- 
lar cylinder. An arbitrary number of blades is 
allowed to be attached to the hub. The blades can 
have arbitrary taper, sweep, dihedral and twist 
angle and can be formed from an arbitrary number of 
different section shapes. The vortex sheet (in the 
case of a circulation that varies along the blade 
span) is assumed to leave the blade from the trail- 
ing edge and continues downstream without allowing 
for the roll-up process. The shape of the vortex 
sheet is arbitrarily prescribed and kept constant 
during. the calculation. Therefore, the vortex 
sheet is allowed to be transparent, that is, it was 
not treated as a stream surface. The iteration 
sweeps start from the line of mesh points connect- 
ing the upstream infinity with the leading edge 
stagnation point (see Fig. 2) and proceed along the 
upper blade surface and then along the lower blade 
surface towards the trailing edge, relaxing <p on 
one line of the mesh points at a time. In this way 
the sweeping direction coincides with the main 
stream direction. Consequently, the artificial 
viscosity introduced by the upstream differencing 
(designated with the superscript H in Eq . (25)) 
will always have positive sign, thus making the 
scheme stable in the regions of locally supersonic 
flow. 


After the iteration converges on the first 
(very coarse) mesh, the values of qp are inter- 
polated onto the next finer mesh (having approxi- 
mately eight times as many points as the first one), 
thus providing an improved initial guess for the 
iterative process on that mesh. The same procedure 


is repeated with the finest mesh after the process 
converges on the second mesh thus resulting in an 
accelerated iterative scheme. 


Preliminary Results 

In order to .test the WIND program for a highly 
compressed relative flow, a two-bladed wind turbine 
rotating at 55 m.p.h. with on oncoming wind speed 
of 18 m.p.h. was studied. The geometric character- 
istics of this rotor are shown in Fig. 6. 

The computation was performed on a single very 
coarse mesh, which consisted of 24x6 mesh cells for 
each two-dimensional plane (Fig. 2). The s panwise 
distribution had 6 mesh cells on the blade surface 
and 2 additional off the blade tip. After 60 itera- 
tions the relative Mach number distribution on tne 
suction side of the blade surface was plotted 
(Fig. 7). This figure indicates that the compress- 
bility effect increases significantly from hub to 
tip with the tip actually operating in transonic 
speed regime. Since such small number of itera- 
tions were performed and the computational mesh was 
so coarse, the results shown in Fig. 7 are prelimi- 
nary ones . 


Summary 

A computer program was developed that numeri- 
cally solves an exact mathematical model for three 
dimensional rotating steady flow through a propel- 
ler-type wind turbine rotor of arbitrary geometry. 
The air is assumed to be mviscid and compressible. 
This work uses the principles of modern computation- 
al aerodynamics and provides designers with a prac- 
tical tool for determining more efficient aero- 
dynamic shapes of wind turbine blades. 
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Figure 3. - Cylindrical computational planes and 
boundary conditions. 


Figure 4. - Geometric transformation sequence. 





Figure 5. - Local isoparametric trilinear mapping. 
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Figure 6. - Spanwise distribution of the blade 
geometric parameters. Airfoil section: 
NACA 2300 series; rotor diameter: 250 ft. 
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